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ABSTRACT 

A 4000 Hz vibration phenomena has been observed during the 
test firings of several space shuttle main engines (SSME) . 
Experimental studies of this phenomena suggest that the problem 
might be associated with vortex shedding from the vanes within 
the LOX tee manifold. The objective of this study was to 
determine the unsteady, feh©ee=Sa?mensionai G ‘ flow associated with 
these vanes by computational methods to identify and better 
understand the 4000 Hz vibration phenomena. A flow solver, FDNS, 
for the turbulent conservation equations was validated for 

predicting high-frequency vortex dynamics and used to predict 2-D 

( r . ® ^ t # 

d-wnens-iona-l and 3“d4me.n>s4©na4^ flows within the LOX tee. 4000 Hz 
excitation oscillations were predicted for some flows and the 
entire 3-dime.aSii«®>na»l flow structure was predicted for LOX tee 
flow. The complexity of the flow was revealed by this analysis, 
and computational methods for predicting these high frequency 
oscillations in future engine systems were established. 
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INTRODUCTION 

A 4000 Hz vibration phenomenon has been observed during the 
test firings of several space shuttle main engines (SSME) . The 
phenomenon is manifested by high acceleration levels in the 
longitudinal direction at the gimbal block location, as shown in 
Fig. 1 for engine 2025 (Ref. 1) . Data from engine accelerometers 
indicate that the vibration is isolated to the gimbal block and 
main injector head. The LOX manifold for the main injector head 
consists of a "hot dog" shaped inlet tee as shown in Fig. 2. The 
purpose of the manifold is to uniformly distribute the flow of 
oxygen. Turning vanes and slots within the LOX tee are used to 
control the flow pattern. Possible fluid dynamic contributions 
to the problem are: turbulence and ballcock instabilities from 

the MOV (these caused major engine problems and required redesign 
of earlier engines) , vorticies caused by the elbow between the 
MOV and the tee inlet, vortex shedding at the edges of the 
turning vanes, and edge tones of the slots. Extensive 
experimental studies have been conducted to understand the 4000 
Hz phenomenon and the flows which contribute to it. This study 
was undertaken to describe, by computational methods, the complex 
flows which occur in this subsystem of the SSME. 

Based on the limited experimental evidence available at the 
time, NASA and Rocketdyne implemented a modification to the LOX 
tee vanes to alleviate the 4000 Hz problem. The opening between 
the two vanes was enlarged, and the trailing edge of the vanes 
was sharpened to shift the frequency of the shed vortex system. 
These modifications reduced the severity of the problem and are 
still in use. However, the value of understanding these flows in 
better detail is evident. 

Detailed review of hot fire test data revealed that the 4000 
Hz oscillation was present in all engines; however, its level 


1 





SECA-TR-90-02 



1 

I 


4120 psla ^ HYDROGEN CAVITY 

196 °R 

69.6 1bm/ft J 

855 Ibm/sec Figure 2. LOX Tee Manifold 

184 ft/sec 



SEC A-TR-90-02 

Was smaller in 

meaSUrS " ents ZLZZZ*".* " bU22er " engines. 

S^rrcjr Mx *- - zrzizz be — 

^orations showed o, d °nie. Th^ i 3v - 

«"9e; the spi ltter va V^**> *>«*■ in the 3600-4000 Hz 
n ° des «e 3500 to 4500 1 T”’ "* t0rsion natural Z 

z::rr ^ d ™ 

f a«. The vibration occurred , t0 hlgher levels on 
^nmgs. since the phenomenon did dUri " 9 M,h power level 

33 '» *» -neitive to TlJT ^ 311 engines, lt 

ll9ht ^n-etric variations. 

small S i°" al CFD a ZsZaZLZt minary inViscid - «*o- 
tl th VOrtices in the „ a * e of °“ Str °" 9 v ° r tical now „ ith 
to the vane locations had a , Vane - “—trie variatT 

sttucture. In f aot , - this vortex ^ 

-no 1 e terally^int o ^the ^ *° 

original position. Th t n °" then ^turning it to it! ^ 
analogous to diffuser „ Z"*"'*"** «» Stable and 
observation has not been determined ' ^ Si3nl ««nce of this 

n zz *» 

This study provides a th»- 
d ynamic analv<?io « three “ d ^ n,e nsiona3 ^ 

toedline from th Zov^T^ “““ 

^itional detailed s^TtT ^ ««=■ 

°f 9 ve° f tW °" dimensio nal blades were mlTZ "*** thS tralli n 9 

£ - :r — 

ss fflechanisns of the 4000 - 

tlne :zz reai 


4 


SECA-TR-90-02 


characteristics of the wake flow. Turbulence models were 
carefully selected, such that, realistic vortex structures at the 
trailing edge of the turning vanes were simulated. 

The computational methodology required to analyze the 
unsteady flow in the LOX manifold tee will be discussed in order 
to establish the critical features which must be correctly 
simulated and the computation techniques which must be used to 
yield good results. Appropriate verifications and meaningful 
sample problem demonstrations were identified and investigated. 
Details of the computer code required and the results of the 
computational investigation are presented in the remainder of 
this report. 

APPROACH 

SECA's study was directed solely toward investigating the 
flowfield features which can excite the 4000 Hz vibration. 
Probably the phenomenon is a fluid/ structures interaction, but a 
sufficiently detailed analysis of the flowfield to identify the 
flow driving the vibration must treat the structure as rigid in 
order to be tractable with state-of-the-art CFD methodology. To 
determine an optimum solution of the problem, the source of the 
oscillation must be determined. Following such a determination, 
a rational decision of whether to change the structural dynamics 
or the flowfield characteristics of the LOX manifold could be 
made. 


The objective of this study was to determine the unsteady, 
three-dimensional, turbulent flow in the LOX manifold tee. 

Ideally, the MOV, the high pressure duct from the MOV to the 
LOX tee, the LOX tee proper, and at least a portion of the LOX 
manifold should be included in the computational simulation. 
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Treating all of these components simultaneously was found to be 
impractical, therefore the three-dimensional cases analyzed were 
started at the inlet to the elbow located just upstream of the 
LOX tee and terminated at the slots between the LOX tee and the 
LOX manifold. In fact, several very complex problems required 
solution before meaningful simulations to this restricted region 
of the flowfield were achieved. Restricting the flow region 
under consideration, as just stated, should not greatly affect 
the engineering utility of the analysis. 

The flowfield analysis presented herein required that 
turbulence effects introduced by vorticial flows in the elbow and 
those arising from interaction with vanes, be evaluated with 
respect to their possible contribution to the 4000 Hz phenomenon. 
If minor geometric variations between engines is indeed the 
explanation of why "buzzers” occur, vane geometry and wall 
boundary conditions within the LOX tee are crucial to the 
analysis. 

To provide an analysis which describes the required flow 
detail, a three-dimensional, unsteady, incompressible, turbulent 
flowfield must be predicted. The solution must be time-accurate, 
which means that pressure waves must be allowed to travel in the 
analysis at the rapid sound speeds typical of the supercritical 
thermodynamic conditions prevalent in the LOX tee. A 
determination of a suitable flow solver was made, and its 
validity was verified. In addition to time accurate pressure 
predictions, boundary condition specifications must not be 
allowed to introduce computational noise into the solution. 

Waves adjacent to boundaries must not be produced by numerical 
error in the algorithm. Proper specification of subsonic 
boundary conditions was made so that such problems were not 
encountered. 
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The critical flow feature, which turbulence models and the 
solution methodology must correctly predict, is the vortex- 
induced vibrations from the trailing edge of the turning vane. 
This places a severe restraint on the type of turbulence model 
and on the grid density required. 

The required simulation was performed by making an 
evaluation of the CFD code to be used, evaluating the effect of 
minor vane geometry variations with a two-dimensional 
calculation, and finally calculating the entire 3-dimensional LOX 
tee flowfield. These steps will be described in the next 
sections of this report. 

FLOW SOLVER SELECTION 


Three concepts are available for obtaining time-accurate 
solutions of the conservation equations for liquid flow. (1) A 
thermal equation of state may be used to introduce the actual 
wave speed for pressure pulses into the system of governing 
equations. (2) An incompressible, constant density analysis may 
be made whereby pressure wave speeds are treated as traveling at 
infinite speeds, i.e. there is no damping due to even small 
density changes. (3) Subiterations between time steps may be 
used to allow artificial compressibility methods to relax to 
obtain time accuracy. The approach proposed for this study was 
to develop a time-accurate CFD algorithm based on a realistic 
equation of state and then apply the method to the LOX tee 
investigation. This development resulted in generating the LWIND 
code, which is presented in Appendix A. However, before the 
validation of this code was completed, two other codes became 
available: INSUP2D and FDNS . INSUP2D was a modification to the 

NASA/Rocketdyne INS code for incompressible flow using artificial 
compressibility, which used subiterations to obtain time accuracy 
(Ref. 3). The equation solved for pressure contains unsteady 
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terms and a parameter to define the artificial compressibility. 
Solution of this equation is hyperbolic in nature and results in 
pressure waves progressing through the computational domain 
during the course of the solution. The speed of these waves 
depend on the value of artificial compressibility. The FDNS code 
(Appendix B) utilizes a pressure based soltuion which solves an 
approximate Poisson equation for pressure. This equation 
contains no unsteady terms, and its solution is elliptic like, in 
that pressure waves have traveled at infinite speed throughout 
the computution domain to satisfy the pressure boundary 
conditions at all points on the boundary of the domain. Rather 
than belabor the value and degree of approximation used, INSUP2D 
and FDNS were both used to solve wake flow problems. If such 
problems could be accurately solved with either or both codes, 
development would be suspended on the LWIND and the more 
developed codes would be used. Before this comparison is 
presented, a discussion of the real fluid thermal equation for 
LOX will be made. 

• The LWIND Code 

To describe unsteady, 3-dimensional, turbulent flow of LOX, 
the interaction of pressure, density and temperature must be 
accurately simulated. Utilizing the parameters: /3 , the 

coefficient of volume expansion, and k, the isothermal 
compressibility, the thermal equation of state, as derived in 
Appendix A, is: 


P = (/3 T + In p - C p )/k 

For reference conditions of pressure equal to 4000 psia, 
temperature of 200 °R, and density of 69.134 65 lb^ft 3 , the 
constants: C p, k, and /? are: 
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C p = 4.59672 
k = 1 . 3948E-5 (psi" 1 ) 

/? = 2.0823E-3 (R' 1 ) 

The values of these parameters were obtained from Ref. 4. In 
general, k and {3 are not constant, but are functions of 
temperature and density. However, the simplified equation of 
state is useful for determining an upwind algorithm for liquid 
flow which is approximately isobaric and isothermal. 

The more direct approach is to determine an accurate 
equation of state, such an equation was developed in Ref. 5, 
based on the methods of Ref. 6. This equation is referred to as 
the HBMS equation of state, it is lengthly, but explicit for 
pressure as a function of density and temperature. The HBMS 
equation is quite accurate for a wide range of liquid and vapor 
flow conditions, including the supercritical. Although the HBMS 
equation gives accurate property values, it is too complex to 
incorporate into an eigenvalue solution for use as a switch in an 
upwind CFD algorithm. The k,(3 equation is adequate for the 
latter purpose. For problems with large density and temperature 
variations, the HBMS equation can be used to modify local « and / 3 
values during the course of the computation. It should be noted 
that artificial compressibility is a numerical parameter which is 
unrelated to fluid properties through an equation of state. This 
means that calculated pressure waves would propagate through the 
medium at a finite, but arbitrary velocity, if artificial 
compressibility methods were used. 

• Evaluation of the INSUP2D and FDNS Codes for Wake Analyses 

The incompressible, time-accurate CFD model, INSUP2D, was 
investigated by comparing computed results to measured vorticial 
flows, specifically, cylinders in cross-flow. Study of the time- 


9 



SECA-TR-90-02 


accurate, high order (5th, 3rd or 1st) upwind INSUP2D code was 
initiated with the intention of using and modifying this code for 
extensive use. A test run for a steady laminar flow past a 
cylinder at. Re = 40 was successfully run with the artificial 
compressibility set equal to 5 and at = AT = 0.25. The 
artificial compressibility is a constant parameter in the 
equation which relates pressure to the divergence of the velocity 
field. The residual (divmax) dropped 4 orders of magnitude in 
1000 steps. Flow at this Re is steady, as was the calculated 
result calculated with the unsteady code. The code changes 
required to generalize the grid to other than cylindrical 
geometries and to add a turbulence model were identified and 
implemented. When a case is run, about half a dozen parameters 
must be specified. As the case is being run, intermediate checks 
are made to insure that these parameters were set at the proper 
levels . 

As stated previously, the laminar wake behind a circular 
cylinder for Reynolds number 40 was predicted to be asymmetric by 
the INSUP2D code. The built-in initial boundary disturbance in 
subroutine BCIMP, used to initiate vortex shedding, rather than 
the grid used, caused the asymmetry. The same case was rerun 
with the initial disturbance turned off. The same steady-state 
input parameters were used. A steady-state solution was obtained 
in 2000 time steps. CPU time was 16 min, 16 sec. The results 
show a much more symmetric wake pattern. See Figs. 3 and 4 for 
stream function and pressure coefficient contour plots, 
respectively. 

Experimental studies of the Re = 40 case have been presented 
by Coutanceau and Bouard (Ref. 7). The parameters reported are 
defined in Fig. 5 and summarized in Table 1. The experiment did 
not conform exactly to infinite flow around a cylinder because 
the flow was contained in a large pipe, therefore the 
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Figure 3. Streamlines at Re = 40 for a Cylinder in Cross-Flow. 
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Figure 4. Pressure Coefficient Contours at Re = 40 for a 
Cylinder in Cross-Flow. 
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Figure 5. Geometrical parameters of the closed wake for a 
Cylinder in Cross-Flow. 


13 


SECA-TR-90-02 


Table 1. Re = 40 Comparisons 


Parameter 

Experiment 

Prediction 

D 

1.0 

1.0 

L 

2.13 

2.2 

a 

0.73 

0.77 

b 

0.59 

0.57 

l 

1.08 

1.08 

■'■max 

X lmax 

0.99 

1.0 

r s 

53.5° 

52.9° 

C D 

1.55* 

1.48 


* from Braza, et al 
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investigators extrapolated their results to eliminate this 
effect. Their extrapolated results were compared to several 
other sets of data from the literature and found to be in good 
agreement with most of the other investigations. Table 1 also 
shows the results of the INSUP2D prediction; the prediction is 
excellent. Experimental results at Re = 10, 15, 20 & 30 were 
also reported by Coutanceau & Bouard. Qualitatively, all of 
these tests produced the same closed wakes with standing 
vorticies. All of the parameters could not be measured for the 
lower Reynolds number cases. For Reynolds numbers lower than 1, 
the standing vorticies would not be expected, but such flows are 
not considered to add to our understanding of the turbulent flows 
in the LOX tee. Additional cases run in the range of the 
referenced data should give comparisons very much like those 
shown in Table 1. 

Cases of unsteady wake flow behind the same cylinder for Re 
= 200 were also simulated using the INSUP2D code. Two cases with 
and without the initial disturbance imposed were studied. The 
same mesh system with grid size of 100 x 120 was used. The 
physical and psuedo-time step sizes assigned were 0.025 and 
1.0E10 respectively. The compressibility parameter was 2500.0, 
and the fifth-order upwind scheme was used. These parameters 
were suggested by Rogers and Kwak. A nearly periodic solution 
was obtained after 2000 time steps (CPU time = 5 hrs, 10 min, 10 
sec) when the initial disturbance was activated. This results in 
a Strouhal number of 0.1942. Figs. 6 and 7 show the final stream 
function plot and pressure coefficient contours. For the second 
case without using the initial disturbance to generate flow 
disturbances, no vortex shedding pattern was predicted after 2000 
time steps. The two cases, with and without initial boundary 
disturbances, were calculated until periodic solutions were 
obtained. With the initial disturbances, the first case produced 
a periodic solution with a Strouhal number of 0.194 in 2000 time 
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Figure 6. Streamlines at Re = 200 for a Cylinder in Cross-Flow. 
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Figure 7. Pressure Cofficient Contour at Re = 200 for a Cylinder 
in Cross-Flow. 
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steps. A time step of 0.025 was required and about 5 hours of 
CPU time on the CRAY-XMP were used. For the second case, the 
flow started to oscillate at about 2000 time steps; it took 
another 2000 steps to reach a perodic solution. The overall CPU 
time for the second case was 7 hrs, 39 min, 46 sec, i.e., 6.89 
sec/time step. It was also found that an average of 20 
subiterations were required per time step. The same Strouhal 
number was obtained for both cases. Rogers & Kwak report a 
Strouhal number of 0.185 for these conditions; they reference a 
value of 0.194 calculated by Lecointe & Piquet. Experimental 
data from Schlichting indicate a value of 0.19. The unsteady 
laminar flow cases were satisfactorily predicted at Re = 200. 

Reynolds number 100, 200, & 1000 cases were studied 
computationally with a laminar model by Braza, Chassaing, & Minh 
(Ref. 8). The reported Strouhal numbers at Re = 200 was 0.2; 
Roshko's experimental value was 0.194. SECA calculated a value 
of 0.194 with INSUP2D. The referenced Re = 1000 case indicated 
that three types of eddies contributed to the wake: (1) those 

associated with the separation point, (2) those associated with 
the free shear layers which separate the free stream from the 
wake, and (3) the classical Karman vorticies. Laminar, 2-D 
calculations are only approximate in this flow regime because the 
wake flow is turbulent and 3-D flow features are of importance. 
Turbulence models can account for both of these effects, if they 
are designed to do so. The turbulent energy spectra in the wake 
of a cylinder at Re = 540 are reviewed by Hinze (Ref. 9) . Braza, 
et al, acknowledge the limitations of the analysis which they 
present, and they point out that the relative importance of the 
eddy structures change as the Reynolds number increases. It is 
known that as the Reynolds number increases, a stable turbulent 
wake is established with a laminar boundary layer ahead of the 
separation point on the cylinder, and at further increases the 
cylinder boundary layer becomes turbulent and the wake decreases 
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in size (Ref. 10) . Classically, the Strouhal number is based on 
the observed vortex structure in the cylinder wake. Between Re = 
5xl0 5 and 4xl0 6 , vortex structure is not experimentally observed 
in the wake of cylinders (Ref. 11) . All of these flow features 
are quite interesting but they can not be addressed 
computationally with a laminar flow analysis. The Re = 200 case 
run with INSUP2D was as successful as could be expected without 
using turbulence models. Other cases for Reynolds numbers up to 
200 can be investigated with the original version of the INSUP2D 
code, which was for laminar, 2-dimensional only. Both of these 
restrictions must be removed before INSUP2D could be used to 
describe the LOX tee flow. 

• Addition of Turbulence Models to the INSUP2D Code 

Algebraic and low-Reynolds number, k-e turbulence models 
were added to INSUP2D to determine if this code could be useful 
for simulating LOX tee flow. 

The Baldwin-Lomax turbulence model (Ref. 12 AIAA Paper 78- 
257) was incorporated into INSUP2D . The cylinder in cross-flow 
was run at a Reynolds number of 10 6 with this turbulence model. 
Streamlines and pressure coefficient contours are shown in Figs. 

8 & 9. The periodicity of the solution is demonstrated with Fig. 
10 which shows the oscillation of the lift and drag coefficients. 
The Strouhal number predicted for this case was 0.29, not the 
0.22 expected; therefore, the Baldwin-Lomax model alone does not 
predict accurate vortex shedding for cylinders. The run time for 
this case was 14.5 hours. 

To improve the turbulence models in INSUP2D, a major effort 
was devoted to incorporating low-Reynolds number turbulence 
models into the INSUP2D code. The turbulence modeling constants, 
source terms, and wall damping functions of Chien's low Reynolds 
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Figure 8. Streamlines at Re = 10 6 for a cylinder in Cross-Flow. 
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Figure 9. Pressure Coefficient Contours at Re = 10 6 for a 
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Figure 10. Lift and Drag Coefficients at Re = 10 6 for a Cylinder 
in Cross-Flow. 
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number two-equation turbulence model were added to the code (Ref. 
13). A turbulence model control parameter, ITURB, was added for 
the selection of turbulence models. The convection terms and 
diffusion terms were also added to the code. Standard upwind 
methods were employed to treat the convection fluxes, and a 
central differencing scheme was used to approximate the diffusion 
terms. These left-hand-side terms were assembled and solved in a 
subroutine called STEPKE. 

The basic features and performance of the eight most well 
known low-Reynolds k-e turbulence models have been reviewed and 
tested by Patel, et al. (Ref. 14). These turbulence models were 
compared for several turbulent boundary layer flows with and 
without pressure gradient effects. All of the two-equation 
turbulence models considered in this study are summarized in 
Table 2. No one model was consistently valid for all the cases 
tested. Patel, et. al., concluded that, in general, the Chien 
and Lam-Bremhorst models performed the best. Low-Reynolds number 
turbulence models have not been developed to the point where they 
are engineering tools, rather they are computational models which 
must be carefully assessed before each application. For example. 
So, et al (Ref. 15) have recently suggested an important 
improvement in treating flows with pressure gradients which is 
applicable to this study. 

The general form of these models can be expressed by: 
k t + uk,, + vky - (y k k x ) x - (i/ k ky) y = P r - e - D Eq.l 

e t + ue x + ve y - (u e e x ) x - {u c e y ) y = (C^P,. - C 2 f 2 e) (e/k) - E 

Eq. 2 
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Table 2. Turbulence Models 

The general form of two-equation k-e turbulence models can be 
written as: 

p D t k = V[ (M+Mt/^)vk] + p (P r -e+D) 

P D t e = V[(M+Mt /o e )V£] + p[(e/k) (C 1 f 1 P t -C 2 f 2 e)+E)] 

where the turbulent viscosity and kinetic energy production rate 
are given by: 

H t =pC (J f^k 2 /e 

P r = (Mt/p ) [ 2 (u x 2 +v y 2 +w z 2 ) + (v x +Uy) 2 + (w y +v 2 ) 2 + (u z +w x ) 2 - (2/3) (u x +v y +w z ) 2 ] 
D and E are defined by the following parameters: 

R t = turbulent Reynolds number = p k z //j.e 

Ry = dimensionless distance from the wall = py/(k)/fi 

y + = another dimensionless distance from the wall 
= py7(TW/p)/M 

Also note: A«7(10j/k/e) 

P + =(M/p 2 0 ( 3 sP)w 
U r =7(r w /p) 

The turbulence model constants, C M , C x and C 2 , are tuned 
against basic turbulent flows (e.g. homogeneous turbulence, wall 
equilibrium conditions, planar and circular jets, etc.)* The 
turbulence Schmidt numbers, a k and a e , are determined based on the 
spreading rate of k and e which satisfy the consistency condition. 

The model constants and damping functions used in several 
k-e models are summarized in subsequent pages of this table. 
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Table 2. Turbulence Models, cont. p.l 


Model 

standard 

extended 

Launder- 



Chen-Kim 

Sharma 


B.C. wall functions 

wall functions 

K = £« = o 

C, 

0.09 

0.09 

0.09 

Cl 

1.44 

1 . 15+0 . 25Min{ 3 , 

P r /e } 1.44 

C 2 

1.92 

1.90 

1.92 

°k 

1.0 

0.8927 

1.0 

O' 

1.3 

1.15 

1.3 


1.0 

1.0 

exp{ -3 . 4 ( 1 +R t /50) " 2 } 

fl 

1.0 

1.0 

1.0 

f 2 

1.0 

1.0 

1-0 . 3exp{ -R t 2 } 

D 

0.0 

0.0 

-2i/ (SyVk) 2 

E 

0.0 

0.0 

2^^ t (3yyU) Z 

Model 

Hassid- 

Hoffman 

Dutoya- 


Poreh 


Michard 


B.C. k„ = e w = 0 

K = e w = 0 

K = = 0 

Cm 

0.09 

0.09 

0.09 

Ci 

1.45 

1.81 

1.35 

C 2 

2.0 

2.0 

2.0 

0 k 

1.0 

2.0 

0.9 


1.3 

3.0 

0.95 

fM 

1-exp {-0. 0015R t } 

exp{ -1 .75/ (1 +R t /50) } 

1-0 . 86exp{ - ( R t / 6 0 0 ) 2 } 

fl 

1.0 

1.0 

l-0.04exp{-(R x /50) 2 }+(A/2y) 

f 2 

l-0.3exp{-R T 2 } 

l-0.3exp{-R x 2 } 

1-0. 3exp{-(R T /50) 2 }-0. 08 (A/y) 

D 

-2i/k/y 2 

-(i//y)9 y k 

-2v(d y Jk) z 

E 

~2v { dyj e ) 2 

0.0 

-C 2 f 2 (eD/k) 2 
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Table 2. Turbulence Models, cont. p.2 


Model : 

Reynolds 

Chien 

Lam- 

Bremhorst 

Kr e w B.c. 

kv,=0, € w =udyyk 

^ W 0 

kv,=0, e„=^3 yy k 


0.084 

0.09 

0.09 

Ci 

1.0 

1.35 

1.44 

c 2 

1.83 

1.80 

1.92 

K 

1.69 

1.0 

1.0 


1.3 

1.3 

1.3 

K 

1-exp { -0 . 0198Ry} 

1-exp { -0 . 0115y + } [1- 

exp{-0. 016 5Ry) ] z (l+2 0. 5/R T ) 

fi 

1.0 

1.0 

1+ ( 0 . 05/ f^) 3 

fz 

fJl-O.Sexpt-RxVs)] 

1-0.22 exp { -R T 2 / 36} 

1-exp { -R T 2 } 

D 

0.0 

-2t/k/y 2 

0.0 

E 

0.0 

—2i> (e/y 2 ) exp{ -0 . 5y + } 

0.0 

Model 

Nagano- 

Lai-So- 



Hishida 

Hwang 


B.C. 

O 

II 

S 

VO 

II 

A? 

K = = 0 



0.09 

0.09 


Ci 

1.45 

1.35 


C 2 

1.90 

1.80 



1.0 

1.0 


<>< 

1.3 

1.3 


■K [i- 

exp { -y + / 2 6 . 5 } ] 2 

1-exp {-0. 0113y + (l-4 . 

■ 372p + ) } 

fi 

1.0 

1.0 


fz 

l-0.3exp{-R T 2 } 

1” ( 2/9 ) exp { -R t 2 / 3 6} 


D 

-2v( dyjk) 2 

-2uk/y 2 


E 

(SyyU) 2 

-2i / (e/y 2 ) exp{-0. 5y + ) 
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where 


P r = i/ t [(Uy + v x ) 2 + 2(u x 2 + v y 2 ) ] , in 2 dimensions. 

"k = v + v J a k 
v t = v + vJo t 
j/ t = C M f„k 2 /e 

D, E, f : , f 2 and f„ are near-wall or low Re damping terms. 

C^, °k > a €> c if and c 2 are turbulence modeling constants. 

There are basically two ways of achieving near-wall damping 
among the turbulence models reviewed in Ref. 14. One group of 
models utilizes the D and E terms for near-wall damping (e.g. the 
Lauder-Sharma and Chien models) . Others employ f x and f 2 terms 
in Eq. 2 to achieve the same effect (e.g. the Reynolds and Lam- 
Bremhorst (LB) models) . The Chien and Lam-Bremhorst turbulence 
models were chosen for further study to represent the two 
different methods of near-wall damping. 

The source terms and near-wall damping terms of the Chien 
and LB turbulence models were coded in subroutine VIST of 
INSUP2D. A new subroutine called STEPKE was created to assemble 
the matrix coefficients of the turbulence equations. In STEPKE, 
the right-hand-side terms of Eq. 1 and 2 are first computed. The 
convection terms were discretized with a first-order upwind 
approximation. A second-order central scheme was employed for 
all diffusion terms. Finally, the free-stream and solid-wall 
boundary conditions for a flat-plate boundary layer flow and a 
circular cylinder wake flow problem were coded. To solve the 2 x 
2 matrix equations, two options, were provided: the original LU 

line-relaxation method and a pointwise over-relaxation matrix 
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method. 

To determine if the numerical formulation and coding were 
correct, a flat plate turbulent boundary layer was analyzed. A 
measured turbulent boundary layer profile (Ref. 9) was imposed at 
the upstream boundary. The computational domain extends 4 
initial boundary layer thicknesses in the cross-stream direction 
and 40 boundary layer thicknesses in the streamwise direction. 
Results using both the Chien and Lam-Bremhorst models shown in 
Figs. 11 and 12 indicate similar boundary layer development. 

With a grid size of 50 x 120 both models gave converged steady- 
state solutions in 1000 time steps. These results also show 2.31 
times increase in CPU time when the turbulence models are 
activated. As indicated by the results of these studies, the 
incorporation of low-Re models into INSUP2D code was successful. 


In order to compute boundary layer type flow problems, the 
INSUP2D code was rearranged to provide the correct geometric 
configuration. The conditions for boundary layer type flows are 
activated by setting the periodic boundary condition flag to 
false in the input data file. For the flat plate mesh system, 
more grid points were clustered near the solid wall, which is 
located at j =1, to provide enough grid points (around 10 
points) inside the viscous sublayer (y + ,10). A subroutine for 
generating the initial turbulent boundary layer profiles (at 
k=kmax) was also created. It was found that the low-Reynolds 
number turbulence models were very sensitive to the inlet 
turbulent boundary conditions since the wall damping terms are 
very sensitive functions of the turbulence quantities. As a 
result, the profiles of the inlet turbulence quantities 
(especially the non-measurable turbulent energy dissipation rate, 
e) must be carefully selected or provided by experimental data. 
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U-VEL0C1TIES 



TURBULENCE-ENERGY 



Figure 11. Boundary Layer Flow with the Chien Low-Re k€ Model at 
X = 40 S 0 . 
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After the two low-Reynolds number turbulence models were 
added to the INSUP2D code and the turbulence models were verified 
to predict boundary layer flow on a flat plate at Re = 10E6, 
flow about a cylinder was analyzed. A solution for pressure 
contours and streamlines is shown in Fig. 13. Predicted 
turbulent kinetic energy is shown in Fig. 14. The free-stream 
was specified to have 1% TKE; the wake flow was found to have 
more than 30%. This value is not realistic, which means that 
this turbulence model requires additional tuning. Because this 
computation was very slow in terms of machine time used, 
developing a three-dimensional solution for the LOX tee manifold 
starting with INSUP2D is not recommended and was not further 
developed. 

SECA's FDNS code has been used to simulate flow about a 
cylinder at Re = 40 and 200, Refs. 16 and 17, respectively. The 
results of these predictions are of comparable accuracy to those 
just presented for the INSUP2D study, but the calculation was 
made in less than 1/10 the time. Therefore, FDNS was selected 
for further development to investigate wake flow behind vanes and 
3-dimensional LOX tee flow. 

BLADE WAKE STUDIES 


To expedite the study of vane wakes required for 
investigating the 4000 Hz problem, detailed two-dimensional 
analyses of flows around blades typical of those in the LOX tee 
were made. However, the time accuracy of such solutions was 
first validated by investigating two vane geometries for which 
experimental data were available. These validation cases were 
studied to determine if the effects of implicit damping in the 
algorithms and/or in the turbulence model overdamped the 
formation of vorticies in the wake of the blades. If the 
numerical flow solvers correctly predicted the vortex shedding 
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Pressure Contours 



Streamlines 


Figure 13. Turbulent Flow for a Cylinder in Cross-Flow with 
Chien's Low-Re Turbulence Model. 
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frequency for these cases, more complex conditions could then be 
considered. 

• Vane Validation Studies 

Predictions of vortex shedding characteristics of rounded 
and sharp trailing edge blades were made. A second-order time- 
centered finite difference Navier-Stokes flow solver (FDNS2D) was 
used in this study. The flowfields around blade configurations 
which were experimentally investigated by A. Haji-Haidari and C. 
R. Smith (Ref. 18), were simulated. In these experiments, the 
water flows, upstream of the trailing edge, are turbulent with a 
Reynolds number of 32,000 (based on blade thickness and an 
average velocity) . Mean velocity and turbulent intensity 
profiles inside the wake were measured using hot film anemometry. 
Flow visualization and recording techniques utilized a hydrogen 
bubble generator and a high-speed video system to determine the 
vortex shedding frequencies. 

Two cases with flow Reynolds numbers of 1,000 and 32,000 
were numerically simulated for each blade in the present study. 
For the low Reynolds number cases, the flow was assumed to be 
laminar. For the high Reynolds number cases, an extended two- 
equation turbulence model with wall function boundary conditions 
was employed for transient turbulent flow computations. A grid 
size of 101 x 73 was used. All of the 2-dimensional vane 
calculations were started from a computational plane normal to 
and cutting through the vane. This procedure greatly reduced the 
number of grid points required for a simulation. The experiments 
were performed in a duct such that bounding walls inclosed the 
vane tested. These walls were used with slip boundary conditions 
to bound the computation domain. In order to provide good 
temporal resolution, a dimensionless time step size of 0.025 
(based on blade thickness and mean velocity) was used for each 


32 


SECA-TR-90-02 


case. Time-accurate periodic solutions were obtained in 5000 
time steps for the Re = 1000 laminar flow past the round trailing 
edge blade starting from an initially uniform flow field. Fig. 

15 (a) -(c) illustrates the solutions of velocity fields, pressure 
contours and stream function contours, respectively, at the end 
of 5000 time steps. An additional 3000 time steps were used to 
obtain a periodic turbulent flow solution starting from the 
laminar flow solution. Fig. 16 (a) -(c) shows the flowfield 
solutions for the turbulent flow case. Results of these two 
cases are quite similar in terms of the vortex shedding patterns 
except that the strength of the vortices are stronger for the 
laminar flow case. Both cases produce Strouhal numbers around 
0.21. The experimental measurements for the turbulent flow case 
agree with this value, indicating that the present flow solver 
was time accurate. 

For the sharp trailing edge cases, however, no periodic 
vortex shedding patterns were predicted even at the end of 6000 
time steps. Pressure oscillations near the tip of the trailing 
edge indicated a very high frequency mode. This suggests that 
much better grid resolution and much smaller time-step size are 
required in order to accurately predict the vortex shedding at 
such high frequencies. This high frequency feature for this 
blade shape was also observed in the experiments (Strouhal number 
= 2.16 - 6.03). Fig. 17 (a) -(c) shows the present flow field 
solutions at the end of 6000 time steps. 

The test case with the cylindrical trailing edge 
geometrically resembles the LOX tee vane trailing edge. This 
test case sheds vortices at about 1 Hz? for LOX manifold flow 
conditions and assuming the same Strouhal number, the frequency 
would be about 1100 Hz. Higher Reynolds numbers, three- 
dimensional flow effects, curvature of the vane, velocity 
imbalance on either side of the vane, and/or an inadequate 
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(b) Pressure Profiles. 
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Figure 15. Laminar Flow over a Vane. 
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■ Figure 15. Laminar Flow over a Vane. 
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(c) Streamlines. 


I 

I 

I 

I 


Figure 17. Turbulent Flow Over a Vane with a Sharp Trailing 
Edge. 
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turbulence model might contribute to the indication of too slow 
an excitation frequency being caused by the flow. 

• LOX Vane Parametric Studies 

Since the blade wake studies were conducted before results 
of the 3-dimensional flow calculations were available, several 
blade configurations with various flow conditions were studied 
parametrically. These simulations were used to determine if 
oscillating vorticies which could cause a 4000 Hz excitation 
could be generated by non-uniform velocity fields and blade 
curvature. 

Two blade configurations were studied at a Reynolds number 
of 10 6 which is representative of that found in the LOX tee. One 
blade was straight; the other had the curvature of the blades in 
the LOX tee. Both of these blades were computed with bounding 
walls which roughly simulated those in the LOX tee. The FDNS2D 
code was used for these calculations with a k-e turbulence model. 
Results for the flat blade case are shown in Fig. 18. The 
predicted Strouhal number for this case was 0.3. The increased 
Strouhal number was a result of the higher Reynolds number and of 
the bounding side walls being closer to the blade than was the 
case for the experimentally tested blade previously described. 
Note this result also differs from that experimentally observed 
for a cylinder for which a Strouhal number does not exist at this 
Reynolds number. Three reasons are suggested for a possible 
explanation of this difference: (1) the turbulence model used 

might not be appropriate, (2) experimentally a time average flow 
is observed, if an ensemble average were observed the vortex 
structure might be present for the flow over the sphere, and/or 
(3) the flow for a blade might be different than that for a 
cylinder. Regardless of the reason (s) for the discrepancy, the 
flows predicted with the k-e model are yielding interesting 
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(a) Pressure and Velocity Oscillations. 



Instantaneous Streamline Locations 


(b) Streamlines. 


Figure 18. Straight Blade Uniform Inlet Velocity (Re 
K-e Model 2000 Time Steps) . 
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results and are as realistic as can be expected without a more 
extensive turbulence data base. For this blade thickness (0.22 
in) and a representative velocity (188 fps) , this corresponds to 
a frequency of 3.07 KHz. Such a frequency is characteristic of 
the experimentally observed wake structure and the oscillating 
lift coefficient, C L . The drag coefficient, C D ; velocity; and 
pressure (in the near wake) oscillate at twice this frequency. 
This doubling of frequencies was also observed by Braza, et al 
(Ref. 8). 

The case with the curved blade is shown in Figs. 19 & 20. 

The frequencies are essentially the same, but the amplitude of 
the pressure and velocity oscillations essentially double. The 
drag coefficient curve shown in Fig. 20 does not clearly indicate 
a doubling of frequency; rather the curve takes a more complex 
shape. Since the lift and drag coefficients are integrated 
values over the blade surfaces, the blade tips may tend to 
oscillate at a higher frequency than the entire blade. In the 
LOX tee, different velocities are expected on either side of the 
blade. Even with uniform velocity on either side of the vanes, 
the pressures and velocities shown in Fig. 18 indicate 
oscillations of about 3.08 kHz, which leads one to speculate that 
more realistic flowfield predictions could increase these values 
to the full 4 KHz in the primary mode without considering 
harmonics. This strongly suggests that the 4 kHz problem might 
be directly simulated with a CFD analysis. 

Also notice that for the curved blade case both the blade 
and the bounding walls are curved and parallel to each other, 
much like the flow in the LOX tee. The experiment and analysis 
reported by Ref. 19 to analyze blade flow in the SSME LOX tee was 
for a curved blade in a straight walled duct. Undoubtedly, the 
duct walls caused the separation which was observed in both the 
experiment and the simulation. SECA does not believe that such 
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Figure 20. Drag coefficients for the curved blade (Re = 
10 6 2000 Time Steps) . 
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flow is representative of that found in the LOX tee. 

To make the 2-dimensional simulations more realistic, cases 
were studied which had different velocities on either side of the 
vane. Specifically, (1) flow over straight and curved blades 
with mismatched velocities on either side of the blade and (2) 
flow over the modified blade configuration of the SSME were 
simulated. 

A flat blade wake case for which the velocities on either 
side of the blade vary by 20 percent, i.e. 207 fps on one side 
and 169 fps on the other, was run at a Reynolds number of a 
million with the FDNS code and a k-e turbulence model. Figure 21 
shows the pressure and velocity oscillations and an instantaneous 
streamline plot. The apparent Strouhal number was 0.27; 
pressures and velocities oscillated at twice that frequency. The 
straight blade with uniform velocity at Re = 10E6 produced lift 
coefficient oscillations of 3078 Hz. The straight blade with 
velocity mismatch of 20% produced lift coefficient oscillations 
of 2828 Hz. Complex modes of the oscillation patterns were 
predicted for this case. 

The curved blade flow with a 20% velocity mismatch was run 
for a Re = 10 6 and the results are shown in Fig. 22. A 
preliminary run was made at a Re = 10 3 to start the turbulent 
flow calculation. The flow is quite different from any blade 
flow that we have previously described. The velocity and 
pressure oscillations appear to be of the same frequency as the 
lift and drag coefficient oscillations; the near wake flow does 
not resemble other flows investigated. The lateral wake 
oscillations appear quite small in amplitude. The wake 
oscillations occur at a frequency of about 3913 Hz. Note: This 

appears to be the source of the 4000 Hz oscillation. 
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The blades in the LOX tee manifold have been modified in 
order to eliminate the 4000 Hz vibration problem in the SSME. 

The trailing edges of the blades have been sharpened and the 
opening between the two blades has been enlarged. The effect of 
the sharpened trailing edge was studied with a two-dimensional 
CFD analysis using FDNS2D. The grid shown in Fig. 23 was used 
for this analysis. Modified vane cases are summarized in Fig. 

24. Uniform flow on either side of the blade was simulated. A 
Re = 10 3 case was run until a periodic solution was obtained, 
then the Re was increased to 10 6 . The periodic solution is shown 
in Fig. 25 for pressure contours and flow streamlines. Notice 
the smoothness of the streamlines and the pressure pattern does 
not indicate vortical flow; the blade shape change is performing 
as expected. Pressure and velocity were monitored at a point 
close to the trailing edge, point A, and were found to almost 
stabilize and then begin to oscillate at a high frequency. 
Although the C L and C D stabilized at about 20 units of 
dimensionless time, about 60 units of time were required before 
pressure and velocity became periodic. This behavior is shown in 
Figs. 26 and 27. Since the amplitude of the C L and C D 
oscillations were small, the pressure and velocity fluctuations 
were believed to be highly localized. Plots of pressure and 
velocity at points B, C, and D, located as shown in Fig. 27, are 
shown in Figs. 28 - 30. These figures indicate that these 
oscillations do die out very near the trailing edge. 

The new trailing edge configuration was further studied by 
analyzing a velocity mismatch case. Again the turbulent flow 
solution was started from a periodic average Re = 10 3 solution. 
The pressure contours and streamlines are shown in Fig. 31 for an 
average Re = 10 6 . Some vortex structure is shown very near the 
trailing edge for this case. The same temporal behavior of 
pressure, velocity, C L , and C D as was found for the uniform flow 
case is seen in the predictions shown in Figs. 32 and 33. 
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(b) Drag and Lift Coefficients 
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Figure 26. Solution History at 
(Re = 10 6 ) . 
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Figure 28. Solution History for a Modified Vane, at Point B. 
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Figure 29. Solution History at Point C for a Modified Vane. 
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Figure 30. Solution History at Point D for a Modified Vane. 


55 





(a) Pressure Contours. 
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t Flow 


(Re = 10 6 ) . 
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(b) Drag and Lift Coefficients. 

Figure 33. Solution History at Point A for Non-Uniform Turbulent 
Flow. 
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The results of the new blade analyses compared to those for 
the old round blade indicate that almost all of the vortical flow 
in the blade wake should be removed. 

The significant 2-dimensional blade wake flow analyses are 
summarized in Table 3. These predictions indicate that 4000 Hz 
oscillations could be predicted in the LOX manifold with this 
computational model. Performing a three-dimensional analysis of 
the LOX manifold using the FDNS3D code was the next step 
performed in this study. 

LOX TEE MANIFOLD ANALYSIS 


• LOX Tee Geometry 

The geometry of the LOX tee manifold from the main LOX valve 
to the slots downstream of the vanes was described with 
analytical expressions in preparation for creating a 3-D 
computational grid for analyzing the flow in the manifold. All 
of the required geometric features were quantitatively 
established from system drawings and from discussions with MSFC 
and Rocketdyne personnel. The analytical expressions were 
incorporated into a computer code to generate the grid for the 
LOX tee. The code was compiled, executed, and verified on a 
microcomputer using standard FORTRAN. 

Figs. 34 through 39 present various views of the geometry as 
represented by a bare minimum 5940 nodes. The actual geometry 
requires between 100,000 and 150,000 nodes. Fig. 34 shows the 
90-degree elbow upstream of the tee. Fig. 35 shows the plane of 
symmetry (Z=0 plane) and illustrates relative postions of the 
elbow, the "ring" connecting the elbow and the tee, and the tee 
itself. Fig. 36 shows two views of the tee with the splitter 
vane. The 0.25 inch radius fillet around the top and bottom of 
the vane is exaggerated due to the coarse grid. Fig. 37 shows 
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Tapered Curved Vane 
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Figure 37. Exit Plane of the Tee. 
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Figure 39. The Hot Dog - Outside View. 



65 



SECA-TR-90-02 


the grid distribution across the exit plane of the tee. Figs. 38 
and 39 show views of the tee and hotdog. For clarity, the 
surface containing the slots has been removed and is shown 
separately. The slots have been darkened in to highlight their 
locations . 

The code was transported to MSFC's CRAY computer where the 
final grid was generated and verified. Half of the symmetric 
geometry was modeled, with 99660 nodes. This distribution 
represents a minimum number of nodes necessary to approximate the 
Loxtee geometry. Stretching of the nodal distributions in all 
three directions was employed in an effort to minimize the number 
of nodes required. 


• 3-D LOX Tee Flow 

In order to reduce the required memory size, for a problem 
requiring a large grid size, FDNS3D was reorganized. Presently, 
FDNS3D can be run with 3.9 mega-words of memory on the Cray-XMP. 
Since some of the metric coefficients must be calculated 
repeatedly, this modification caused a two percent increase in 
CPU time for every time step. However, due to the nature of the 
non-iterative time marching scheme employed in the FDNS3D code, 
reasonably efficient turn-around time was still realized. 

Initially, computation of a 3-D laminar flow field for the 
Loxtee was performed. The flow Reynolds number, based on the 
inlet mean velocity and the vane thickness, was assumed to be 
1000. Results of this laminar flow problem were used as an 
initial flowfield for a turbulent flow case. The code ran at a 
speed of 1.75E-04 sec/grid/time-step on the Cray-XMP 
supercomputer which means that about 17 sec/time-step were 
required for the grid size used. A nearly periodic flowfield 
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solution, based on the monitoring values of velocity and pressure 
downstream on the vane trailing edge, was obtained in 1250 time 
steps. The preliminary results of this laminar flow study are 
shown in Figs. 40 and 41 and indicate a pressure oscillation 
frequency close to 4400 Hz. Figure 40 illustrates the 
perspective views of an instantaneous flow pattern inside the LOX 
tee geometry. It can be seen from Fig. 40(c) that a large 
velocity difference on both sides of the vane was predicted. 

This may have contributed to the high frequency pressure 
oscillation (compared with the case of uniform flow past a vane) . 

The velocity profiles upstream of the vane trailing edge, 
which were shown previously in the 2-D blade study to have 
substantial effects on the vortex shedding frequency in the wake 
of the vane, are shown in Fig. 42. Much higher flow speed was 
predicted along the concave side of the vane than that on the 
convex side of the vane. This effect is shown in Fig. 42 for the 
velocity profiles at three vane thicknesses upstream of the 
trailing edge. These velocity differences on the both sides of 
the vane can create large shear stress effects on the vortex 
shedding frequency in the wake. Results of this 3-D analysis can 
provide useful inflow boundary conditions for a more complete 2-D 
vortex shedding analyses. 

The computation of a turbulent flow with Reynolds number 
about one million, inside the 3-D LOX tee geometry was 
accomplished. An extended two-equation turbulence model was used 
in the analysis. Several modifications to the outlet and wall 
boundary conditions were made in order to eliminate their 
influences on the final vortex shedding predictions. First, the 
outlet mass flow conservation condition was modified to let the 
flow speed corrections be proportional to the flow rate in each 
slot. The wall function formulation was also modified to reduce 
sensitivity to the near-wall grid skewness. However, all these 
modifications were not sufficient to stabilize the solution for 
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Figure 40. Perspective View for the Instantaneous Flow Pattern 
in the LOX Tee. 
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(C) 


Figure 40. Perspective View for the Instantaneous Flow Pattern 
in the LOX Tee, Continued. 
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Figure 42. U-Velocity Upstream of Vane Trailing Edge. 
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this turbulent flow case. Therefore, the mesh near the cap end 
was modified to reduce the grid skewness near the corner points 
in the computational domain. 

A new mesh system with grid size: 157 x 19 x 33, was 

generated by modifying the cap end of the tee to reduce the 
degree of grid skewness near that region, when this modification 
was made, the turbulent flow case became very stable. Since only 
a small amount of flow passed through the region at the end of 
the tee, this modification should not have significent influence 
on the wake flow solution downstream of the vane. An extended 
two-equation turbulence model was also employed in this 
investigation. 

The same initial guess of the flowfield used for the laminar 
flow case was utilized to start the turbulent flow computation. 
The inlet turbulence intensity and turbulence length scale were 
assumed to be 5 percent of inlet flow velocity and 3 percent of 
the inlet duct diameter, respectively. It took 1600 time steps 
before a nearly periodic solution was obtained. A non- 
dimensional time step size (based on the inlet velocity and the 
vane thickness) of 0.02 was used in this case. 

Figs. 43(a) — (d) illustrate four flow pattern projections 

for J=4, 8, 11 and 16 planes respectively. Large recirculation 
zones on the upper surface of the tee and on the convex surface 
of the vane were predicted. Except for the stronger 
recirculation, the basic features of the flow were very similar 
to the laminar flow solutions. The flow has gone through a sharp 
turn before exiting from the downstream slots. The entire flow 
structure was highly three-dimensional, as can be seen from six 
cross-sectional views (1=61, 71, 81, 91, 93 and 101 planes) of 
secondary flow patterns given in Figs. 44(a) - (f). Two eddies 
with counter-clockwise rotation. Fig. 44(a), were formed upstream 
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(e) 



(f) 



1=101 


Figure 44. Turbulent LOX Tee Flow in The Cross-Plane, Continued. 
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of the vane trailing edge (the vane trailing edge was located at 
1=67) . A strong shear layer was then generated downstream of the 
vane trailing edge, Fig. 44(b). This pattern persists until a 
prevailing flow through the exit slots is established, Figs. 

44(c) - (f). 

The computation of a turbulent flow inside a half-domain LOX 
tee geometry indicated that a periodic solution was obtained in 
2400 time steps. The periodicity was based on the time history 
of the integrated lift forces on the LOX tee vane. This result 
showed a Strouhal number of 0.25, or an oscillation frequency of 
2485 Hz. The predicted low oscillation frequency may be 
attributed to the basic assumption of the present simulation that 
a symmetry boundary condition was imposed along the center plane 
between two vanes. This may have omitted one important feature 
of the LOX tee flow; namely that a bi-stable flowfield 
downstream of the vanes may exist. Experimental observations of 
the water flow tests suggest that the 4000 Hz oscillations may 
have originated from the vane to which the flow was attached 
after a bi-stable flow pattern was established. 

Three extensions of the half-plane simulation were made. 
First, a complete LOX tee geometry was incorporated in the 
computation. This increases the grid size to around 200,000 
nodes, in order to run this problem within 4M words memory size, 
a multi-zone algorithm was implemented by using the SSD device on 
the Cray-XMP. Two blocks of SSD areas were used to store the 
grid and flow variables of two zones. For each time step, data 
for each zone were read from SSD in sequence. Solutions of each 
zone were then stored back to SSD after every time integration. 
This procedure was tested and found to be successful. However, 
due to the necessity of doing SSD input/ output , the CPU time 
required for each time step has increased by a factor of 4 . 
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Three hundred time steps were completed using this method. 

At the end, the flowfield was still symmetric with respect to the 
symmetry plane. Figs. 45(a) - (d) show the velocity fields on 4 
J-planes. The integrated forces on both vanes were also found to 
be coincidental for each time step. This suggested that further 
runs would not be able to produce a bi-stable flowfield. The 
initial guesses of the flowfield and the outlet mass flow rate 
distributions, which were symmetric at the beginning, may have 
prevented the development of a bi-stable solution. 

An alternative to the above approach could also be employed 
to study the effect of the bi-stable flowfield on the vortex 
shedding frequency. An excess mass flow can be assigned along 
the symmetry plane downstream of the vane of the half-domain. It 
can be assumed that the entire mass flow through the gap between 
the vanes is attached to one of the vanes. This approach may 
give reasonable results with more efficient computer turn-around 
time. This approach was not tested, rather a fully coupled, 
entire tee domain simulation was performed on the NASA/Ames Cray. 

In order to have a better turn-around time in the 
computation of a full 3-D turbulent flow around the LOX tee vanes 
using a multizone approach, the NASA/Ames Cray-2 and Cray-YMP 
supercomputers were used to carry out the calculations. The 
Cray-2 and Cray-YMP computers can handle a much larger grid size 
than Marshall's Cray-XMP. A two-block mesh system with a grid 
size of 2 x (157 x 19 x 33) was generated using the grid 
generator. This size of mesh system requires about 8 mega-words 
of memory when using the FDNS code with two-equation turbulence 
models. A new zonal interface routine was incorporated in the 
FDNS code to provide faster data updating along the interface 
without using the solid state device (SSD) . This zonal interface 
interpolation procedure is performed in parallel with the matrix 
solver in the FDNS. This two-zone method was running at speeds 
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of 75 sec/time-step and 25 sec/time-step on the Cray-2 and the 
Cray-YMP computers respectively. These are about 3 and 9 times 
faster than the case using the two-zone SSD on the Cray-XMP. 

Since much better turn-around was experienced using the Cray-YMP 
machine, most of the results presented here were computed using 
the Cray-YMP computer. Eleven consecutive runs (a total of 1699 
time steps) with a dimensionless time step size of 0.02 
(normalized with the vane thickness and the average inlet 
velocity) were performed. Subsequently, an additional run of 517 
time steps with a dimensionless time step size of 0.01 was 
completed . 

With 98,439 nodes representing one-half of the flow domain 
(a zonal interface is located between two vanes) , the grid 
resolution is considered to be coarse compared with that which 
was used for the 2-D cases. It would be difficult to obtain 
accurate time-resolved vortex shedding structures with the 
current grid density. The main objective of this 3-D case is 
therefore aimed at the general flow patterns around the LOX tee 
vanes. 

Figs. 46 shows the final result of the 0.02 time step case 
mentioned above; the calculated lift coefficient and monitoring 
values of pressure and velocity at a point inside the near-wake 

of one of the vanes with respect to the normalized time are shown 
in this figure. It is clear from these results that the 
solutions have not yet reached a quasi-steady-state vortex 
shedding pattern. This may be caused by the coarse grid in the 
wake regions or by disturbances from boundary conditions imposed 
along the exit boundary (which is not far away from the vane 
trailing edges) . However, results of this investigation have 
indicated different lift forces on the two LOX tee vanes which 
also indicate non-symmetric flow patterns on the downstream side 
of the vanes. The same features were not found in the previous 
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Figure 46. Calculated Monitoring Parameters for the Full 3-D LOX Tee Flowfield 
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investigations using two-zone SSD approach on the Cray-XMP. This 
phenomenon is consistent with the experimental observations in a 
water flow test program. Fig. 47 and 48 illustrate several 
cross-sections of the predicted flow pattern at the end of the 
calculation. The predicted velocity distributions upstream of 
the trailing edge of the vane indicate more attached flow along 
the suction surface, in the plane shown on p. 86. Again, large 
velocity differences on both sides of the vane were predicted in 
this investigation. But, the velocity difference magnitudes are 
less than that of the previous one-zone solution with symmetric 
boundary conditions. This is also consistent with the water flow 
test results and the assumption made in the calculations for 2— D 
vane configurations. 

In summary, based on the 2-D results, half a million grid 
points would be required to resolve the vortex shedding structure 
in the LOX tee. Therefore, the 200,000 grid points used for the 
3-D calculation indicate only the general features of the 
flowfield. These general flow features are useful for making a 

thorough 2— D vane wake analysis. Due to the near proximity 
of the surface on which exit conditions are specified to the 
vanes in the tee, exit boundary conditions should be further 
investigated to more accurately simulate the LOX tee flow. 
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CONCLUSIONS AND RECOMMENDATIONS 

Based on the computational results reported herein and the 
apparent success of the proposed fix to the 4000 Hz problem which 
NASA and Rocketdyne are currently using, SECA offers the 
following conclusions and recommendations for this investigation. 

1 • The FDNS CFD code was verified to be useful in analyzing 
unsteady oscillating phenomena at least up to the 6 kHz 
range using wall function boundary conditions and a 
sufficiently dense computational grid. 

2 . The two-dimensional blade wake investigation demonstrated 
that flow disturbances at the 4000 Hz frequency level could 
originate in the wake of the splitter vanes. Furthermore, 
analyses of this type could be used to avoid flow conditions 
and geometric configurations which might cause similar 
problems in future engine designs. 

3. Full 3-dimensional CFD solutions for LOX tee flow were 
generated. These solutions revealed the general features of 
the flowfield, but current computer size and speed 
limitations prevent a sufficiently detailed analysis to 
identify wake fluctuation features. The general flow 
features are suitable for use as boundary conditions for 
detailed 2-dimensional wake analyses. 

4 . Even though the physics simulated in a CFD solution is 
correct, complex 3-dimensional problems, like the LOX tee 
flow, must still be addressed by looking in detail at small, 
select regions of the overall flowfield. Carefully designed 
parametric CFD investigations of unsteady, internal liquid 
flows can yield invaluable design information. 
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APPENDIX A - LWIND 


1. INTRODUCTION 

This document describes LWIND, a code developed by SECA, 
Inc., to provide solutions to complex liquid flow problems. 

LWIND solves the Navier-Stokes equations using an equation of 
state for liquids that accurately defines the liquid in the flow 
regime of interest. The code currently operates in 2-D or 
axisymmetric modes with options for inviscid, laminar or 
turbulent flows. Both Baldwin-Lomax and k-e turbulence models 
have been included in the code. In addition to the algorithm 
discussed below, the code has a grid generator, a preprocessor 
and a graphics package. 

The following sections present the derivations of the liquid 
equation of state and thermodynamic properties and the windward 
algorithm employed by the code. 

2. THERMODYNAMIC PROPERTIES OF LIQUIDS 

This section discusses the derivations of the equation of 
state, enthalpy and internal energy for any liquid. These 
relationships are to be used in a Navier Stokes solver and, 
therefore, must be consistent with the continuity, momentum and 
energy equations. Namely, the density must be the inverse of the 
specific volume, the internal energy and enthalpy must differ 
only by P/p, and the eigenvalues of the set of equations must 
involve only the velocity and speed of sound. For instance, many 
text books assume p is constant in the derivation of the equation 
of state. This leads to an inconsistency between density and 
specific volume which, in turn, results in inconsistencies 
between internal energy and enthalpy and extraneous terms in the 
eigenvalues. The following derivations result in a set of 
relations for the liquid thermodynamic properties which is 
consistent throughout the analysis. 

2.1 Equation of State 

The thermal equation of state for any fluid may be expressed 
in a differential form as 

dp = (0p/ST) p dT + (0p/0P) T dP 

where p is the density, P is the pressure and T is the 
temperature. But, by definition 

( 3p/3T ) p = -p/3 

where 0 is the coefficient of volume expansion and 
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(3p/3P)rp = pK 

where k is the isothermal compressibility. Substituting these 
definitions into the equation of state and rearranging yields 

d p/p = -0dT + /cdP 

It is at this point that the assumption of constant p results in 
the aforementioned inconsistency. Instead, assume that 0 and k 
are constant over the range of interest for a given problem and 
are evaluated at reference conditions p r , P r , and T r . Integrating 
the equation of state yields 

In p = — j9T + kP + C p 

where C p is a constant of integration evaluated at the reference 
conditions. That is 

Cp = In p r - «P r + ,8T r 

For multi-component mixtures of miscible liquids the density 
of the fluid is the sum of the densities of the components. Care 
should be taken in the assumption of miscibility of the 
components. If the components a are mixed at the same pressure 
and temperature, the equation of state becomes 

In p = In ( - TI/3 a + PX« a ) + C p 
a a 

where 0 = I0 a , k = l K a r an< 3 ^p ^ or fixture only. 

a a 


The static pressure is obtained directly from the equation 
of state 

P = (j9T + In p - Cp ) / K 

If the reference conditions are judiciously chosen from 
thermodynamic property data for any liquid, the above equation of 
state will usually be valid over a sufficient range of pressures 
and densities. 

2.2 Enthalpy 

The enthalpy of any fluid may be derived using the same 
process. The differential form of the equation for enthalpy is 
given by (Ref. Sears , T herm o dymanics , p 169) 

dH = (3H/3T)p dT + (3H/3P) t dP 
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= CpdT + [v - T(Bv/ST) p ]dP 

where H is the enthalpy, v is the specific volume, and C p is the 
specific heat at constant pressure. Some mathematical 
manipulation using the Chain Rule leads to 

dH = [C p - ( jS 2 vT//c) ]dT - d v/k + (jS//c)d(vT) 

where, referring to Sears, p. 149, 

C p - ( 0 2 vT/ k ) = C p - (C p - C v ) = c v 

It can be shown (Sears, p. 150) that the specific heat at 
constant volume C v is independent of the specific volume v if the 
pressure is a linear function of the temperature. Referring to 
the equation for pressure above, it can be seen that pressure is 
indeed linear in temperature. Therefore, 

C p - ($ 2 vT/k) s C p - ( ,8 2 v r T / k ) 

Since C p varies little it can be assumed constant. These 
relations reduce the equation for enthalpy to a form which is 
easily integrated. 

dH = [Cp - ( jS 2 v r T//c ) ]dT - dv//c + ( jS / «: ) d(vT) 

Integrating the above equation yields 

H = C p T - (0 2 v r T 2 /2 k) - v/k + j3vT /k + C H 
where Cp is a constant of integration 

C H s H R “ c p T r + ((S 2 T r 2 /2 np r ) + 1/Kp r - ( (ST r /Kp r ) 

Now the equation for enthalpy becomes 
H = C a T - 1/kp + jST /kp + Cp 
where C Q s 0.5(C p + C v ) 

For multiple components the intensive enthalpies are 
additive, that is 

pH = £ p a H a 

a 

The enthalpy equation becomes 

pC a T - 1/k + (ST /k + pC H = I[p a (C a ) a T - 1/K a +Tj3 a /A: a + P a (Cp) a ] 

a 
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Comparison of the two sides of the above relation infers the 
following relations for multiple component liquids. 


P°p = I 


a 


o 

< 

II 

M 

P a ( c v ) 

a 


pC H = l 

p«( c h) 


a 


2.3 Internal Energy 

The equation for internal energy may be derived using the 
same method, starting with the differential form (Sears p. 148) 

dE = ( 3E/ 3T ) v dT + (3E/3v) T dv * C v dT + ( jST /k - P ) dv 

Applying manipulations similar to those for enthalpy and 
integrating yields 

E = CpT -(/S 2 T 2 /2 K p r ) + ( vln v - v + C p v)/* + C £ 

E = C g T - 1/Kp + fiT/np - P/p + C £ 

Comparing this equation to the equation for enthalpy indicates 
that relationship H = E + P/p is satisfied if C £ = C H , 
therefore the equation for internal energy becomes 

E = C a T - 1 / Kp + |8T/*:p - P/p + Cjj 

2.4 Speed of Sound 

A relationship for the speed of sound will be required for 
the eigenvalue analysis to follow. The equation for the speed of 
sound is (Sears p. 155-156) 

C 2 = 1 /pk s 

where k s is the adiabatic compressibility. Denoting the ratio of 
specific heats as y , the adiabatic compressibility becomes 

k s = K/r 

and the speed of sound relationship becomes 
C 2 = r/pic 
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2.5 Derivatives 

The formulation of the Jacobian matrices in the following 
section will require taking the derivatives of pressure and 
enthalpy. These derivatives were derived realizing that the 
partial derivatives of the independent variables with respect to 
each other are zero. 

2.5.1 Pressure 

Starting with the relationship for pressure 
P = (/ST + In p - Cp) / k 

the derivative of pressure with respect to 6 becomes 
SP/Sp = C 2 - a ( Kip - ) 

where a is defined as $/KpC v . Also, 

9P/9pu = <3 ( -u/pC v ) /k - - au 
9P/9 pv = - av 
9P/9 pEp = a 

2.5.2 Total Enthalpy 

SHp/ 9 p = [C 2 - a(Hp - q 2 ) + Hp]/p 
9Hp/9pu = - au /p 
9Hp/9pv = - av/p 
3Hp/9pp = ( a + l)/p 

The relations derived in this section will be used in the 
next section to provide the thermodynamic properties of the 
liquid for the windward algorithm. 


3 . WINDWARD ALGORITHM 

An explicit windward algorithm using a finite volume flux 
splitting technique has been developed for solving the 
conservation equations for a liquid. The explicit method was 
chosen because of the need to solve real time problems. Total 
conservation in the integral sense is guaranteed by employing the 
elemental flux integrals to form the spatial derivatives found in 
the conservation equations. The viscous terms in the flow 
equations are resolved using standard central differences and are 
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not included in the following discussion of the derivation of the 
algorithm. The energy equation is included to permit the 
treatment of heat transfer or liquids with different energy 
levels . 

For 2-D and axisymmetric flow, the conservation equations 
may be cast in the form: 

Q t + F x + Gy + o H/y = 0 


Q = 

[P, 

pu, pv, pE T ] T 

F = 

[pu. 

p U 2 + p, puv, puH T ] 

G = 

[pv, 

puv, pv 2 + P, pvHj] 

H = 

[pV, 

puv, pv 2 , pvH T ] T 


where t is time; x & y are Cartesian coordinates; u & v are 
velocity components in the x & y directions, respectively; p is 
density; P is pressure; E T is total specific internal energy; H 
is total specific enthalpy and <7 = 1 for axisymmetric flow and 
= 0 otherwise. 

The conservation equations may be transformed to local , 
computational coordinates ( £ , n ) where A £ = A 77 = 1 . 

(Y ff /J)Q t + (Y a F)^ + ( y° G)„ + <xH = 0 

where : 

J = ^x V Y - £ y r; x 

*x = j Yt? 

V x = - J Y£ 

£y = ~ J x r} 

Vy = J Xg 

F = [pU, puU + £ x P/J, pvU + iy P/J, pUH T ] T 
G = [pV, puV + >7 X P/J, pvV + Vy P/J, pVH T ] T 
H = [ 0 , 0 , -P/J, 0] T 
U = (e x /J)u + (£ y /J)v 
V = (D x /J)u + (r? y /J)V 

The vectors F and G can be locally linearized by first 
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obtaining the Jacobian matrices using the formulas 


- || Q e = A 

3Q 


G 


n 


SO Q = B Q 
0Q 


The mathematical treatment of the Jacobian matrices A and B 
is identical, therefore by letting D represent A or B, n 
represent f/J or rj/J, and U = n x u + nyV, the D matrix becomes 


0 

n x 

n y 

0 

-uU + n x 0 2 

U + (l-a)un x 

uny - avn x 

an x 

-vU + nyfl 2 

vn x - auny 

U + (l-a)vny 

an y 

U(0 2 - H t ) 

**T n x " 

H T n y - avU 

( 1+a ) U 


where 

6 2 = 3P /dp = C 2 - a(H^. - q 2 ) 
q 2 = u 2 + v 2 


The D matrix can been shown to be equivalent to the 
corresponding matrix for an ideal gas by noting that, for an 
ideal gas 


a = 7 - 1 
IS = 1/T 
K = 1/P 

e 2 = 0.5(7 - i)q 2 . 


The eigenvalues for the matrix D are obtained by setting the 
determinant of D to zero and solving for the roots. The 
eigenvalues are equivalent to those for gases, namely 

X - [U, U, U+nC, U-nC] 
where n = (n x 2 + riy 2 ) 0 ' 5 . 

The eigenvector matrix for D is obtained by solving the 
following relationship for each of the four eigenvalues. 
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where V 

II 

3 

X 

< 

1 

3 

•< 

C 

m l 

= u + n x C/n 

m 2 

= v + nyC/n 

m 3 

= H t + UC/n 

m 4 

= u - n x C/n 

m 5 

= v - HyC/n 

m 6 

= H T - UC/n 


and C is the speed of sound. The determinant of the eigenvector 
martix is 

IM d ! = 2 C 3 n/a 

The inverse of M D is obtained using the standard formula 
m d _ 1 = M ji/ im d' 
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where M^j = (-l) a+ ^c^j and c^j are the minors of M^. Then 

,2 


m d _1 


= («/C 2 ) 


H T -q‘ 

VC 2 /an 2 
( f? 2 -UC/n) /2a 
( 0 2 +UC/n ) /2a 


u 


n y C 2 /an 2 
(n x C/an-u) /2 


-n x C 2 /an 2 
( n y C/ an-v ) /2 


-1 

0 

1/2 


-(n x C/an+u) /2 - ( n y C/an+v ) /2 1/2 

It has been verified that is equal to the identity matrix. 

The matrix D may be spectrally decomposed into components 
corresponding to each of the eigenvalues by using the spectral 
theorem. Since two of the eigenvalues are identical the 
decomposition results in only three components. These components 
are obtained by performing the following matrix multiplications 


D = Mr 


or 


U 


1 


0 


0 

1 

+ (U+nC) 

0 

+ (U-nC) 

0 

0 


1 


0 

0 


0 


1 


D — + D2 + Dg 

Di = ( a U/C 2 ) x 


M, 


-1 


(C 2 -e 2 )/a 

u 

V 

-1 

(n x UC 2 /n 2 -u0 2 ) /a 

u 2 +n y 2 C 2 /an 2 

uv-n x n y C 2 /an 2 

-u 

(n y UC 2 /n 2 -v0 2 ) /a 

uv-n x n y C 2 /an 2 

v 2 +n x 2 C 2 /an 2 

-v 

(U 2 C 2 /n 2 -0 2 H T )/a 

Udi-n y VC 2 /an 2 

v ai+n x VC 2 /an 2 

- w 


where w = q 2 - 0 2 /a 
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D 2 = a ( U+nC ) /2C 2 


W 1 

w 2 

w 3 

1 

W^m^ 

W 2 m^ 

w 3 m l 

m 

W2m 2 

W 2 m 2 

W 3 m 2 

m 

w l m 3 

W 2 m 3 

W 3 m 3 

m 



where Wj = (0 2 - UC/n )/a 
W 2 = (n x C/an - u) 

W 3 = (nyC/an - v) 


D 3 = a ( U-nC ) /2C 2 


w 4 w 5 W 6 1 

W 4 m 4 W 5 m 4 W 6 m 4 m 4 

w 4 m 5 w 5 m 5 w 6 m 5 m 5 

_ W 4 m 6 W 5 m 6 W 6 m 6 m 6 


where 

W 4 = ( 0 2 + UC/n ) /a 
W 5 = -(n x C/an + u) 

W 6 = -(nyC/an + v) 

The quantity DQ n now becomes ( D 2 + D 2 + D 3 )Q n where Dj has eigen 
direction U, D 2 has eigen direction U + nC and D 3 has eigen 
direction U - nC. Recall that DQ n represents either 



It 

<*i>« + 

(f 2 

u + 


or 

G v = 

(Si)„ + 

CNJ 
<C D 

U + 

(G 3 ), 


The components Dj , 

d 2 

and 

D 3 may be simplified by applying 


the chain rule to replace derivatives of (pE T ) with derivatives 
of P. The equation for the differential form of the internal 
energy, dE, from section 2.3 may be expanded to the differential 
of the total energy, dEj, to give 

dE T = C v dT - (/ST /k - P)d p/p 2 +qdq 
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But 


dT = /cdP//S - dp//3 p 


d ( pE-p ) — p d E <ji + E-pdp 
d { pq) = pdq + qdp 


Combining the above four relations yields 

d(pE T ) = dP/a - 0 2 dp/a + ud(pu) + vd(pv) 


Substituting the above equation into the three components of D 
eliminates many of the terms in the matrices. Another 
simplification involves using the chain rule to replace (pu) n and 


(pv) n with (pU) n and (pV) 


1 n 


These two operations reduce the 


complicated component matrices into much more simplified forms, 
namely : 


DiQn 


where Q* = 

[p, pU, pV, P] T 

and 


1 

0 

0 

-1/C 2 

= u 

n x U/n 2 

0 

-n y /n 2 

-u/C 2 


n y u/n 2 

0 

n x /n 2 

-v/C 2 


u 2 /n 2 -e 2 

/a 0 

V/n 2 

-w/C 2 



-U 

1 

0 

n/C 

D* 2 = U+nC 

-Um 2 

m l 

0 

m 1 n/C 

2nC 

-Um2 

m 2 

0 

m 2 n/C 


-Um 3 

m 3 

0 

m 3 n/C 


U 

-1 

0 

n/C 

D* 3 = U-nC 

Um 4 

— m^ 

0 

m 4 n/ 

2nC 

Um 5 

-m 5 

0 

m 5 n/C 


Um 6 

-m e 

0 

m 6 n/C 
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D*jQ* n can be expanded by noting that 
(pU) n = (pU) n 

(puU) n = -uUp n + (u + n x U/n 2 ) (pU) n - (n y U/n 2 ) ( pV) n 
(pvU) n = -vU p n + (v + n y U/n 2 )(pU) n + (n x U/n 2 ) (pV) n 
(pUH T ) n = -U(H t + e 2 /a)p n + (H t + U 2 /n 2 )(pU) n 
+ (UV/n 2 ) ( pV) n + U ( a + l)P n /a 

Substituting the above relations into puts the actual 

fluxes back into the algorithm. The algorithm then becomes the 
fluxes plus biases in the windward directions. D* 1 Q* n becomes 

D *lQ*l = F n +[Up n -(pU) n ] 




" " 

1 


U/n 

a l 

-nP n /C 2 

b 3 

a 2 


b 2 

a 3 

1 

b 3 


where F n is either of the original flux vectors, F or G , and 
= u + n x U/n 2 
a 2 = v + n y U/n 2 
a 3 = H T + U 2 /n 2 
b 2 = (uU + n x C 2 )/n 
b 2 = (vU + n y C 2 )/n 
b 3 = U(H t +C 2 )/n 

Defining several convenient groupings of parameters 


[Up n 


(pU) n ] 







1 


U/n 

a l 

f 2 = nP n /C 2 

b l 

a 2 


b 2 

a 3 


b 3 
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U/n 


1 

f 3 = [U p n - (pU) n ]/C 

b l 

f 4 = nP n /C 

a l 


b 2 


a 2 


b 3 


a 3 


then D* lQ * n = (F 2 ) n = F n + fj - f 2 

D*2Q 3 *'n anc ^ ®*3 ( 5*n can a ^- so be expanded to yield 

D* 2 Q* n = (F 2 ) n = 0 . 5 ( f 2 - f 1 ) + 0 . 5 ( f 4 - f 3 ) 

D* 3 Q* n = ( F 3 ) n = 0.5(f 2 - f 2 ) - 0.5(f 4 - f 3 ) 

Again the eigen directions of Fj , F 2 and F 3 are U, U + nC 
and U - nC, respectively. 

Once the flux components have been evaluated, the windward 
algorithm determines which nodes in the element receive each of 
the flux components. The windward algorithm is stated as: 

F n = e l( F l)n + e 2 ( F 2 ^ n + e 3 ( F 3 )n 
where: e^ = 1 if U > 0 

= 0 if U < 0 

e 2 = 1 if U + nC > 0 

= 0 if U + nC < 0 

e 3 = 1 if U - nC > 0 

0 if U - nC < 0 

For supersonic flow, which is unlikely for liquid flows, e^ = e 2 
= e 3 , so that 

F n = e l< F 1 + F 2 + F 3 ) n 

■ ejt F + fj - f 2 + 0.5 (f 2 - f 2 + f 4 - f 3 ) 

+ 0.5 (f 2 -f 1 - f 4 + f 3 )] 

= e l F = e l (F) n 

For subsonic flow where U < nC the e's may be redefined as: 
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e l 

= 0.5(1 

+ 

£ 1 ) / 

where 

£ 1 

e 2 

= 0.5(1 

4* 

£ 2 ) ' 

where 

£ 2 

e 3 

= 0.5(1 

+ 

£ 3 ) ' 

where 

£ 3 


In subsonic flow any node in the 
e 3 equal to unity and the other 
node 


= SIGN ( U ) 

= SIGN( U+nC ) = SIGN (n) 

= SIGN(U-nC) = SIGN(-n) = - e 2 

flow grid will have either e 2 or 
equal to zero, therefore for any 


F n = 0.5[(l+ £l )F + e 1 (f 1 -f 2 ) + e 2 ( f 4 -f 3 )] 

By simply evaluating the signs on U for supersonic flow and for U 
and n for subsonic flow the algorithm automatically determines 
the proper directions in which to difference the three components 
of the flux vectors. The method by which the differencing is 
accomplished in LWIND is described below. 

LWIND employs a numbering system convention which applies to 
each element in the flow field regardless of orientation. This 
numbering system is detailed below. 


Typical Element 


4 C 3 


D 

V 


B 


1 £ A 2 


Nodes : 1-4 

Faces : A - D 

Elements do not need to be 
rectangles; sketch indicates 
numbering system only. 


The position vector for any point in the element is P = [x, y] T 
or 

P = ( 1-e ) ( 1-77 ) P a + H1-J?)P 2 + + (l-£)r?P 4 

P e = ( l-J? > (P 2 " Pi) + V(P 3 ~ P 4 ) 

P„ = (l-e)(P 4 - P 2 ) + *(P 3 - P 2 ) 

IAi = x £ y v - X 7J y £ = 1/3 


In order to understand the physical interpretation of the 
equations used in the algorithm, first note that the outward- 
pointing normals 
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on each face of 
Face 
A 
B 
C 
D 


N = N x i + N y j 
the element are: 

—x $y 

y 2 -yi x i _x 2 

y3~y2 x 2 _x 3 

y 4 -y 3 x 3 _x 4 

yi-y 4 x 4 - x i 


Evaluation of the derivatives and at each node of the 
element yields 


Node 

x ? 

— 

ye 

X r? 

y^ 

1 

K 

ro 

\ 

K 

!-» 

y 2 -yi 

x 4 " x l 

y 4 -yi 

2 

x 2 -x l 

y 2 -yi 

x 3 -x 2 

Y3 _ y2 

3 

x 3 -x 4 

y3-y 4 

x 3 -x 2 

Y3-Y2 

4 

x 3 -x 4 

y3-y 4 

x 4 - x i 

y 4 -yi 


But from before 


tx = J V v V x = - J Vf 

£y = ~ J Xjj 7 ] y = J X£ 


so that 


Node 

*x 

£y 


Vy 

1 

- J 1 N XD 

~ Jl NyJ) 

- Jl n xA 

- J 1 N yA 

2 

3 2 N xB 

J 2 N yB 

” J 2 n xA 

- J 2 Ny A 

3 

J 3 N xB 

J 3 Nyg 

J 3 N xC 

3 3 ^ yC 

4 

- j 4 n xB 

- J 4 Ny D 


J 4 N yC 
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Therefore v£/J and vtj/J are area normals at the 2 faces 
associated with each node in the positive £ and rj directions. U 
is the dot product of the velocity vector and the area normal 
vector. Consequently, the local computational coordinate form of 
the conservation equations represents the fluxes through the 
faces of the element and are, in fact, identical to the flux 
integral form of the conservation equations. LWIND departs from 
standard finite difference methods at this point by employing the 
integral form of the equations and actually integrating the 
fluxes through each element. This results in a finite volume 
integral approach using a windward algorithm. 

The above table reveals that and £ y involve only the 
normals on faces B and D of the typical element while r? x and r; v 
only involve the normals on faces A and C. Consequently, 1 

F is evaluated only on faces B and D while G is evaluated only on 
faces A and C. 

The adopted convention for the typical element, with its 
numbering system for the nodes and faces relative to the f and rj 
directions, dictates the following for e 2 = SIGN(n) . The outward- 
pointing normal for face B points in the direction of positive £ 
while the outward-pointing normal for face D points in the 
negative £ direction. Therefore, 

£ 2 = 1 for nodes 2 & 3 for 

= -1 for nodes 1 & 4 for F^ 

Likewise, the outward-pointing normal directions for faces A and 
C dictate that 

&2 = 1 for nodes 3 & 4 for 

= -1 for nodes 1 & 2 for G^ 

If U a is defined as the average U in either the | or rj 
direction, recalling that U is the dot product of the velocity 
and area normal vectors, then, since the area normal vectors on 
opposite faces have opposite signs, the following is true: 

e 1 = SIGN ( U a ) for nodes 2 & 3 for F^ 

= -SIGN(U a ) for nodes 1 & 4 for F^ 

e 1 - SIGN ( U a ) for nodes 3 & 4 for G^ 

= -SIGN(U a ) for nodes 1 & 2 for 

Applying the above conventions for and e 2 ' to the equation for 
F n and noting that each node on a face receives half of F n 
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yields : 

For nodes 2 & 3 (for F ^ ) or nodes 3 & 4 (for ) 


or 


F + = 

0 . 2 5 ( (1 + £ 2 )F + £ 1 (f 1 

" f 2 ) + 

(f 4 

- f 3 >] 

For 

the opposite nodes, 




F“ = 

0 . 25 [ (1 - £ : )F - £ 1 (f 1 

- f 2 > - 

(f 4 

- f 3 n 

4. 

F- = 

0 . 25 [ (1 + £ j ) F ± 62(^1 

- f 2> ± 

(U 

- f 3 >] 


= 0.25F n ± 0.25(e 1 F n + B] 


where 


® — £ 1 ( f 2 - f 2 ) + (^4 ~ ^3) 

The flux equation consists of a center differenced term involving 
F^ or G^ only and a bias term B which is windward differenced but 
sums to zero over the element, thus ensuring the integral 
conservation of the flow equations over the element. 

The bias term B can be simplified using the following 
groupings of parameters: 



Then 


B = 


£ 2 C - U/n 
e jU/n - C 

[ Up n - (pU) n ]b 1 /C - nP n b 2 /C 2 
[Up n - (pU) n ]b 2 -nP^/C 

b 3 

b 3 u + b 4 n x /n 
b 3 v + b 4 n y /n 
b 3 H T + b 4 U/n 


The above derivation results in a finite volume, totally 
conservative, windward algorithm consisting of center differenced 
fluxes with a windward bias involving no time-consuming matrix 
multiplications . 
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Appendix B - FDNS 

The Navier Stokes conservation equations to be solved, along 
with k-e turbulence model equations, are given below in 
curvilinear coordinates. 

J'Udpq/dt) = dl-pV.q + MGijOq/3^) J/3^ + S q 

where Ui = 1, u, v, w, h, k, and e, respectively, for mass, 
momentum, enthalpy, turbulent kinetic energy, and turbulent 
kinetic energy dissipation conservation. J, U A , and are given 
by 


J = <pU,rt,n/3 (x,y,z) 

U, = (Uj/J) (a^/axj) 

Gu = (3ei/ax k ) (aej/ax k )/j 

Also, n = + Mt)/ a q i s the effective viscosity when the 

turbulent eddy viscosity is used to model turbulent flows. The 
turbulent eddy viscosity is n t = pCjf/e; C M and o q are turbulent 
modeling constants. Wall functions are used to reduce the number 
of grid points required very near no-slip walls. 

The source terms are given by: 


0 


-Px + V[M(Uj) x ] 


S, 


= J* 1 


-Py + ?[MUj)y] 
-p, + V[M(Uj) r ] 
* 


p(Pr - O 

(p€/k) (C x P r - C 2 e + C 3 P r 2 / e ) 


An upwind scheme is used to approximate the convective terms 
of the momentum, energy, and continuity equations ; the scheme is 
based on second and fourth order central differencing with 
artificial dissipation. First order upwinding is used for the 
turbulence equations since the parameters involved are all 
positive quantities. Different eigenvalues are used for weighing 
the dissipation terms depending on the conserved quantity being 
evaluated, in order to give correct diffusion fluxes near wall 
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boundaries (Ref. Bl) . This procedure is different num 
proposed by other investigators (Refs. B2 - B4) in which the sum 
of the absolute value of the convective velocity and the local 
sound speed is used to weigh the dissipation terms. For 
simplicity, consider fluxes in the ^-direction only. That is: 

8F/d( “ 0.5(F i+1 - Fw) “ (d i+1/2 ~ di- 1 / 2 ) 

A general form of the dissipation term is given as follows. 
dj+ 1/2 = 0*5[€i |pu| 

] i+1/2 (q i+ i ~ qi) + [ £ 2 ( l~£i) MAX (0.5ASp(|u|,|v|), 

2 | pU | } + £ 3 A s ] i+ 1/2 (qi~i ” + 3q i+1 — qi+ 2 ) 

where: dj = REC (upwinding control parameter) and d 3 = 0.005 

Different values for e lf e 2r and e 3 are used for the continuity, 
energy and momentum equations, as shown in Table 1. 


Table 1. Dissipation Parameters 


Momentum & Energy 


Continuity 


e 2 0.015 

€3 0 


0 

0 

d 


3 


To maintain time accuracy, a time-centered, time-marching 
scheme with a multiple pressure corrector algorithm is employed. 
In general, a noniterative time-marching scheme is used for time 
dependent flow computations; however, pressure corrector 
subiterations can be used if necessary. The pressure corrector 
scheme is described as follows. A simplified momentum equation 
was combined with the continuity equation to form a pressure 
correction equation. This equation is: 

dp Ui/at “ - vp' 
or in discrete form: 

Ui ' ~ - p ( At/ p ) Vp ' (1) 

where p represents a pressure relaxation parameter ( p = 10 is 
typical). The velocity field in the continuity equation is then 
perturbed to form a correction equation. That is: 

VfpUi) 1 * 1 = VC^fUi" + u^)] =0 


or, 
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V(pU/) = - V(puj n 

Substituting equation (1) into (2) , the following pressure 
correction equation is obtained. 

- V(0At Vp') = - V ( pUj ) n (3) 

Once the solution of equation (3) is obtained, the velocity field 
and the pressure field are updated through equation (1) and the 
following relation: 

p = p + p' 

To ensure that the updated velocity and pressure fields satisfy 
the continuity equation, the pressure correction procedure is 
repeated several times (usually 4 times is sufficient) before 
marching to the next time step. This constitutes a multi- 
corrector subiteration procedure. 

For the description of cavitating flow, an additional 
conservation equation for the gas phase is required. This 
equation contains a source term, pVl t , which is evaluated as 
follows: Derivation of the bubble growth rate two-phase model 

follows the method given in (Ref. B5) , the basic equations used 
in this model are the Clausius-Clapeyron equation and the 
simplified one-dimensional momentum and energy equations. The 
final expression can be written as: 

dm/dp = (v fg /h fg ) (m - CT/h fg ) 

where m, p, v, h, C and T represent the mass fraction (or 
quality) of the bubbly phase, pressure, specific volume, 
enthalpy, specific heat and fluid temperature respectively. The 
subscript £g denotes differences due to phase change. 

The rate of production of heterogeneous nuclei per unit 
volume as proposed in (Ref. B6) is: 

dn/dt = N n (2a /Kn ) 0,5 exp[-(16jra 3 6/3/cT)/ (p g - p f ) 2 ] 

where: N n , number of heterogeneous nucleation sites per unit 

volume ; 

o, surface tension of oxygen; 

M, mass of a mulecule; 
k , Boltzmann's Constant; 

6, function of bubble contact angle, (assumed to be 
unity in this study) ; 
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p g ; vapor pressure at temperature, T; ano 
p £ ; fluid pressure. 

By combining these equations, an expression for the bubble 
growth rate results: 

dr/dt = [K f R 2 T f Vh f6 2 p 6 2 (3Kt/7O 0 5 ] ( Pg 2 / Pgr ) m(p g /p gr) 
where: K £ , thermal conductivity of fluid; 

K, thermal diffusivity of fluid; 

p gr , pressure inside a bubble of radius, r; 

R, gas constant. 

Finally, the source term for the mass fraction equation can 
be written as: 

pW t = p R (l-a) (4?r/3) (U/as) (dn/dt) (dr/dt) 3 

where s is an average cell size and the volume fraction, a, can 
be calculated from the following homogeneous relation: 

a = m(p £ /p g )/[l + m(p £ /p g - 1)] 
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